Water Observations from Space (WOfS)¶

Using Geoscience Australia WOfS algorithm.

This notebook requires the "WOfS Environment", which needs to be installed on the first run of this notebook by executing the following cell. Once done, set the environment in the top-right corner.

In [1]:
# Install the WOfS environment - Only needed on the very first execution of this notebook
!../tools/install_wofs.sh
WOfS Environment kernel already installed

Table of contents¶

  • Data parameters - Lake Tempe
  • Set proxy parameters to access data locally
  • Import dependencies and initialise datacube
  • Display the region of interest
  • Load and display the DEM
  • Load the data
  • Classify WOfS
  • Water Observations Summaries
    • Wet counts
    • Clear observations counts
    • Wet frequencies

Data parameters - Lake Tempe¶

In [2]:
product = "landsat8_c2l2_sr"
longitude = (119.8242517, 120.0350519)
latitude = (-4.2013799, -3.9445384)
time = ('2020-01-01', '2020-12-31')
output_crs = "EPSG:32650"
resolution = (30, -30)

# Where to save the DEM fetched in ODC
DEM_PATH = "/home/jovyan/dems/srtm_lake_tempe.tif"

Set proxy parameters to access data locally¶

In [3]:
from os import environ

environ["AWS_HTTPS"] = "NO"
environ["GDAL_HTTP_PROXY"] = "easi-caching-proxy.caching-proxy:80"
print(f'Will use caching proxy at: {environ.get("GDAL_HTTP_PROXY")}')
Will use caching proxy at: easi-caching-proxy.caching-proxy:80

Import dependencies and initialise datacube¶

In [4]:
import sys
sys.path.append("../../hub-notebooks/scripts")

from pathlib import Path

import xarray as xr
import rioxarray

from app_utils import display_map
from wofs.virtualproduct import WOfSClassifier
In [5]:
from datacube import Datacube
from datacube.utils.aws import configure_s3_access

configure_s3_access(
    aws_unsigned=False, 
    requester_pays=True, 
);

dc = Datacube()
In [6]:
# Display available products
products_info = dc.list_products()
products_info
Out[6]:
name description license default_crs default_resolution
name
copernicus_dem_30 copernicus_dem_30 Copernicus 30m Digital Elevation Model (GLO-30) None None None
landsat5_c2l2_sr landsat5_c2l2_sr Landsat 5 Collection 2 Level-2 Surface Reflect... None None None
landsat5_c2l2_st landsat5_c2l2_st Landsat 5 Collection 2 Level-2 UTM Surface Tem... CC-BY-4.0 None None
landsat7_c2l2_sr landsat7_c2l2_sr Landsat 7 USGS Collection 2 Surface Reflectanc... None None None
landsat7_c2l2_st landsat7_c2l2_st Landsat 7 Collection 2 Level-2 UTM Surface Tem... CC-BY-4.0 None None
landsat8_c2l1 landsat8_c2l1 Landsat 8 Collection 2 Level-1 (top of atmosph... CC-BY-4.0 None None
landsat8_c2l2_sr landsat8_c2l2_sr Landsat 8 Collection 2 Surface Reflectance, pr... None None None
landsat8_c2l2_st landsat8_c2l2_st Landsat 8 Collection 2 Level-2 UTM Surface Tem... CC-BY-4.0 None None
landsat9_c2l1 landsat9_c2l1 Landsat 9 Collection 2 Level-1 (top of atmosph... CC-BY-4.0 None None
landsat9_c2l2_sr landsat9_c2l2_sr Landsat 9 Collection 2 Surface Reflectance, pr... CC-BY-4.0 None None
landsat9_c2l2_st landsat9_c2l2_st Landsat 9 Collection 2 Level-2 UTM Surface Tem... CC-BY-4.0 None None
lpdaac_nasadem lpdaac_nasadem NASADEM Merged DEM Global 1 arc second V001 None None None
nasa_aqua_l2_oc nasa_aqua_l2_oc NASA MODIS-Aqua L2 Ocean Color, regridded to W... None EPSG:4326 (0.01, 0.01)
nasa_aqua_l2_sst nasa_aqua_l2_sst NASA MODIS-Aqua L2 Sea Surface Temperature, re... None EPSG:4326 (0.01, 0.01)
s2_l2a s2_l2a Sentinel-2a and Sentinel-2b imagery, processed... None None None

Display the region of interest¶

In [7]:
display_map(x=longitude, y=latitude)
Out[7]:
Make this Notebook Trusted to load map: File -> Trust Notebook

Load and display the DEM¶

We load the SRTM DEM data from ODC using the lpdaac driver. The Lake Tempe area is mostly flat so the DEM can be made optional.

In [8]:
# Ignore warnings in output
import warnings
from sqlalchemy.exc import SAWarning
warnings.filterwarnings("ignore", category=FutureWarning) 
warnings.filterwarnings("ignore", category=SAWarning) 

from os import environ
from cartopy.crs import PlateCarree
from datashader import reductions
import holoviews as hv
import hvplot.xarray
import matplotlib.pyplot as plt

dem = dc.load(
    product="lpdaac_nasadem", 
    latitude=latitude,
    longitude=longitude,
    output_crs="epsg:4326", 
    resolution=(-1/3600, 1/3600),
)
elevation = dem.elevation
In [9]:
options = {
    'title': 'Elevation',
    'width': 800,
    'height': 500,
    'aspect': 'equal',
    'cmap': plt.cm.terrain,
    'clim': (0, elevation.max().values.item()),    # Limit the color range depending on the layer_name
    'colorbar': True,
    'tools': ['hover'],
}
plot_crs = PlateCarree()
elevation.hvplot.image(
     x = 'longitude', y = 'latitude',         # Dataset x,y dimension names 
     crs = plot_crs,
     rasterize = True,                        # If False, data will not be reduced. This is slow to load but all data is loaded.
     aggregator = reductions.mean(),          # Datashader calculates the mean value for reductions (also first, min, max, las, std, mode)
     precompute = True,                       # Datashader precomputes what it can
    ).opts(**options).hist(bin_range = options['clim'])
Out[9]:
In [10]:
dem_path = Path(DEM_PATH)
dem_path.parent.mkdir(parents=True, exist_ok=True)
elevation.rio.to_raster(dem_path)

Load the data¶

Load the data from ODC and rename bands as needed by the WOfS classifier.

In [11]:
# Ignore SAWarning in output
measurements = ['blue', 'green', 'red', 'nir', 'swir1', 'swir2', 'pixel_qa']
data = dc.load(
    product=product,
    longitude=longitude,
    latitude=latitude,
    time=time,
    output_crs=output_crs,
    resolution=resolution,
    measurements=measurements,
    dask_chunks={'time': 1},
)

data = data.rename({
    "blue": "nbart_blue",
    "green": "nbart_green",
    "red": "nbart_red",
    "nir": "nbart_nir",
    "swir1": "nbart_swir_1",
    "swir2": "nbart_swir_2",
    "pixel_qa": "fmask",
})
data
Out[11]:
<xarray.Dataset>
Dimensions:       (time: 19, y: 951, x: 785)
Coordinates:
  * time          (time) datetime64[ns] 2020-01-09T02:10:26.574371 ... 2020-1...
  * y             (y) float64 -4.65e+05 -4.65e+05 ... -4.366e+05 -4.365e+05
  * x             (x) float64 8.371e+05 8.37e+05 ... 8.136e+05 8.136e+05
    spatial_ref   int32 32650
Data variables:
    nbart_blue    (time, y, x) uint16 dask.array<chunksize=(1, 951, 785), meta=np.ndarray>
    nbart_green   (time, y, x) uint16 dask.array<chunksize=(1, 951, 785), meta=np.ndarray>
    nbart_red     (time, y, x) uint16 dask.array<chunksize=(1, 951, 785), meta=np.ndarray>
    nbart_nir     (time, y, x) uint16 dask.array<chunksize=(1, 951, 785), meta=np.ndarray>
    nbart_swir_1  (time, y, x) uint16 dask.array<chunksize=(1, 951, 785), meta=np.ndarray>
    nbart_swir_2  (time, y, x) uint16 dask.array<chunksize=(1, 951, 785), meta=np.ndarray>
    fmask         (time, y, x) uint16 dask.array<chunksize=(1, 951, 785), meta=np.ndarray>
Attributes:
    crs:           EPSG:32650
    grid_mapping:  spatial_ref
xarray.Dataset
    • time: 19
    • y: 951
    • x: 785
    • time
      (time)
      datetime64[ns]
      2020-01-09T02:10:26.574371 ... 2...
      units :
      seconds since 1970-01-01 00:00:00
      array(['2020-01-09T02:10:26.574371000', '2020-01-25T02:10:22.588145000',
             '2020-02-10T02:10:17.788921000', '2020-02-26T02:10:13.879802000',
             '2020-03-13T02:10:07.386121000', '2020-03-29T02:09:58.450496000',
             '2020-04-14T02:09:51.264053000', '2020-04-30T02:09:43.329680000',
             '2020-05-16T02:09:40.092585000', '2020-06-01T02:09:44.725513000',
             '2020-06-17T02:09:54.647535000', '2020-07-03T02:10:02.482792000',
             '2020-08-20T02:10:18.209098000', '2020-09-05T02:10:25.859941000',
             '2020-09-21T02:10:31.337497000', '2020-10-23T02:10:34.865932000',
             '2020-11-24T02:10:34.546460000', '2020-12-10T02:10:36.006848000',
             '2020-12-26T02:10:33.500205000'], dtype='datetime64[ns]')
    • y
      (y)
      float64
      -4.65e+05 -4.65e+05 ... -4.365e+05
      units :
      metre
      resolution :
      30.0
      crs :
      EPSG:32650
      array([-465045., -465015., -464985., ..., -436605., -436575., -436545.])
    • x
      (x)
      float64
      8.371e+05 8.37e+05 ... 8.136e+05
      units :
      metre
      resolution :
      -30.0
      crs :
      EPSG:32650
      array([837075., 837045., 837015., ..., 813615., 813585., 813555.])
    • spatial_ref
      ()
      int32
      32650
      spatial_ref :
      PROJCS["WGS 84 / UTM zone 50N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32650"]]
      grid_mapping_name :
      transverse_mercator
      array(32650, dtype=int32)
    • nbart_blue
      (time, y, x)
      uint16
      dask.array<chunksize=(1, 951, 785), meta=np.ndarray>
      units :
      reflectance
      nodata :
      0
      crs :
      EPSG:32650
      grid_mapping :
      spatial_ref
      Array Chunk
      Bytes 27.05 MiB 1.42 MiB
      Shape (19, 951, 785) (1, 951, 785)
      Count 38 Tasks 19 Chunks
      Type uint16 numpy.ndarray
      785 951 19
    • nbart_green
      (time, y, x)
      uint16
      dask.array<chunksize=(1, 951, 785), meta=np.ndarray>
      units :
      reflectance
      nodata :
      0
      crs :
      EPSG:32650
      grid_mapping :
      spatial_ref
      Array Chunk
      Bytes 27.05 MiB 1.42 MiB
      Shape (19, 951, 785) (1, 951, 785)
      Count 38 Tasks 19 Chunks
      Type uint16 numpy.ndarray
      785 951 19
    • nbart_red
      (time, y, x)
      uint16
      dask.array<chunksize=(1, 951, 785), meta=np.ndarray>
      units :
      reflectance
      nodata :
      0
      crs :
      EPSG:32650
      grid_mapping :
      spatial_ref
      Array Chunk
      Bytes 27.05 MiB 1.42 MiB
      Shape (19, 951, 785) (1, 951, 785)
      Count 38 Tasks 19 Chunks
      Type uint16 numpy.ndarray
      785 951 19
    • nbart_nir
      (time, y, x)
      uint16
      dask.array<chunksize=(1, 951, 785), meta=np.ndarray>
      units :
      reflectance
      nodata :
      0
      crs :
      EPSG:32650
      grid_mapping :
      spatial_ref
      Array Chunk
      Bytes 27.05 MiB 1.42 MiB
      Shape (19, 951, 785) (1, 951, 785)
      Count 38 Tasks 19 Chunks
      Type uint16 numpy.ndarray
      785 951 19
    • nbart_swir_1
      (time, y, x)
      uint16
      dask.array<chunksize=(1, 951, 785), meta=np.ndarray>
      units :
      reflectance
      nodata :
      0
      crs :
      EPSG:32650
      grid_mapping :
      spatial_ref
      Array Chunk
      Bytes 27.05 MiB 1.42 MiB
      Shape (19, 951, 785) (1, 951, 785)
      Count 38 Tasks 19 Chunks
      Type uint16 numpy.ndarray
      785 951 19
    • nbart_swir_2
      (time, y, x)
      uint16
      dask.array<chunksize=(1, 951, 785), meta=np.ndarray>
      units :
      reflectance
      nodata :
      0
      crs :
      EPSG:32650
      grid_mapping :
      spatial_ref
      Array Chunk
      Bytes 27.05 MiB 1.42 MiB
      Shape (19, 951, 785) (1, 951, 785)
      Count 38 Tasks 19 Chunks
      Type uint16 numpy.ndarray
      785 951 19
    • fmask
      (time, y, x)
      uint16
      dask.array<chunksize=(1, 951, 785), meta=np.ndarray>
      units :
      bit_index
      nodata :
      1
      flags_definition :
      {'snow': {'bits': 5, 'values': {'0': 'not_high_confidence', '1': 'high_confidence'}}, 'clear': {'bits': 6, 'values': {'0': 'not_clear', '1': 'clear'}}, 'cloud': {'bits': 3, 'values': {'0': 'not_high_confidence', '1': 'high_confidence'}}, 'water': {'bits': 7, 'values': {'0': 'land_or_cloud', '1': 'water'}}, 'cirrus': {'bits': 2, 'values': {'0': 'not_high_confidence', '1': 'high_confidence'}}, 'nodata': {'bits': 0, 'values': {'0': False, '1': True}}, 'qa_pixel': {'bits': [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15], 'values': {'1': 'Fill', '2': 'Dilated Cloud', '4': 'Cirrus', '8': 'Cloud', '16': 'Cloud Shadow', '32': 'Snow', '64': 'Clear', '128': 'Water', '256': 'Cloud Confidence low bit', '512': 'Cloud Confidence high bit', '1024': 'Cloud Shadow Confidence low bit', '2048': 'Cloud Shadow Confidence high bit', '4096': 'Snow Ice Confidence low bit', '8192': 'Snow Ice Confidence high bit', '16384': 'Cirrus Confidence low bit', '32768': 'Cirrus Confidence high bit'}, 'description': 'Level 2 pixel quality'}, 'cloud_shadow': {'bits': 4, 'values': {'0': 'not_high_confidence', '1': 'high_confidence'}}, 'dilated_cloud': {'bits': 1, 'values': {'0': 'not_dilated', '1': 'dilated'}}, 'cloud_confidence': {'bits': [8, 9], 'values': {'0': 'none', '1': 'low', '2': 'medium', '3': 'high'}}, 'cirrus_confidence': {'bits': [14, 15], 'values': {'0': 'none', '1': 'low', '2': 'reserved', '3': 'high'}}, 'snow_ice_confidence': {'bits': [12, 13], 'values': {'0': 'none', '1': 'low', '2': 'reserved', '3': 'high'}}, 'cloud_shadow_confidence': {'bits': [10, 11], 'values': {'0': 'none', '1': 'low', '2': 'reserved', '3': 'high'}}}
      crs :
      EPSG:32650
      grid_mapping :
      spatial_ref
      Array Chunk
      Bytes 27.05 MiB 1.42 MiB
      Shape (19, 951, 785) (1, 951, 785)
      Count 38 Tasks 19 Chunks
      Type uint16 numpy.ndarray
      785 951 19
  • crs :
    EPSG:32650
    grid_mapping :
    spatial_ref

Classify WOfS¶

In [12]:
transform = WOfSClassifier(c2_scaling=True, dsm_path=DEM_PATH)
# Compute the WOFS layer
wofl = transform.compute(data)
wofl
Out[12]:
<xarray.Dataset>
Dimensions:      (y: 951, x: 785, time: 19)
Coordinates:
  * y            (y) float64 -4.65e+05 -4.65e+05 ... -4.366e+05 -4.365e+05
  * x            (x) float64 8.371e+05 8.37e+05 8.37e+05 ... 8.136e+05 8.136e+05
  * time         (time) datetime64[ns] 2020-01-09T02:10:26.574371 ... 2020-12...
    spatial_ref  int32 32650
Data variables:
    water        (time, y, x) uint8 dask.array<chunksize=(1, 951, 785), meta=np.ndarray>
Attributes:
    crs:      EPSG:32650
xarray.Dataset
    • y: 951
    • x: 785
    • time: 19
    • y
      (y)
      float64
      -4.65e+05 -4.65e+05 ... -4.365e+05
      units :
      metre
      resolution :
      30.0
      crs :
      EPSG:32650
      array([-465045., -465015., -464985., ..., -436605., -436575., -436545.])
    • x
      (x)
      float64
      8.371e+05 8.37e+05 ... 8.136e+05
      units :
      metre
      resolution :
      -30.0
      crs :
      EPSG:32650
      array([837075., 837045., 837015., ..., 813615., 813585., 813555.])
    • time
      (time)
      datetime64[ns]
      2020-01-09T02:10:26.574371 ... 2...
      units :
      seconds since 1970-01-01 00:00:00
      array(['2020-01-09T02:10:26.574371000', '2020-01-25T02:10:22.588145000',
             '2020-02-10T02:10:17.788921000', '2020-02-26T02:10:13.879802000',
             '2020-03-13T02:10:07.386121000', '2020-03-29T02:09:58.450496000',
             '2020-04-14T02:09:51.264053000', '2020-04-30T02:09:43.329680000',
             '2020-05-16T02:09:40.092585000', '2020-06-01T02:09:44.725513000',
             '2020-06-17T02:09:54.647535000', '2020-07-03T02:10:02.482792000',
             '2020-08-20T02:10:18.209098000', '2020-09-05T02:10:25.859941000',
             '2020-09-21T02:10:31.337497000', '2020-10-23T02:10:34.865932000',
             '2020-11-24T02:10:34.546460000', '2020-12-10T02:10:36.006848000',
             '2020-12-26T02:10:33.500205000'], dtype='datetime64[ns]')
    • spatial_ref
      ()
      int32
      32650
      spatial_ref :
      PROJCS["WGS 84 / UTM zone 50N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32650"]]
      grid_mapping_name :
      transverse_mercator
      array(32650, dtype=int32)
    • water
      (time, y, x)
      uint8
      dask.array<chunksize=(1, 951, 785), meta=np.ndarray>
      Array Chunk
      Bytes 13.53 MiB 729.04 kiB
      Shape (19, 951, 785) (1, 951, 785)
      Count 2469 Tasks 19 Chunks
      Type uint8 numpy.ndarray
      785 951 19
  • crs :
    EPSG:32650
In [13]:
# Uncomment the following line to display WOfS for each time slice. This may take a few minutes
# wofl.water.plot(col="time", col_wrap=5);

Water observations summaries¶

The Water Observations Summaries based on https://github.com/opendatacube/odc-stats/blob/develop/odc/stats/plugins/wofs.py are made up of:

  • count_clear: a count of every time a pixel was observed (not obscured by terrain or clouds)
  • count_wet: a count of every time a pixel was observed and wet
  • frequency: what fraction of time (wet/clear) was the pixel wet
In [14]:
# Rename dimensions as required
wofl = wofl.rename({"x": "longitude", "y": "latitude"})
In [15]:
from odc.algo import safe_div, apply_numexpr, keep_good_only

wofl["bad"] = (wofl.water & 0b0111_1110) > 0
wofl["some"] = apply_numexpr("((water<<30)>>30)==0", wofl, name="some")
wofl["dry"] = wofl.water == 0
wofl["wet"] = wofl.water == 128
wofl = wofl.drop_vars("water")
for dv in wofl.data_vars.values():
    dv.attrs.pop("nodata", None)
In [16]:
# Ignore warnings triggered by time slices without data at all
warnings.filterwarnings("ignore", message="divide by zero encountered in true_divide") 
warnings.filterwarnings("ignore", message="invalid value encountered in true_divide") 

wofl.wet.plot(col="time", col_wrap=5);
In [17]:
# Helper frunction from https://github.com/opendatacube/odc-stats/blob/develop/odc/stats/plugins/wofs.py
def reduce(xx: xr.Dataset) -> xr.Dataset:
    nodata = -999
    count_some = xx.some.sum(axis=0, dtype="int16")
    count_wet = xx.wet.sum(axis=0, dtype="int16")
    count_dry = xx.dry.sum(axis=0, dtype="int16")
    count_clear = count_wet + count_dry
    frequency = safe_div(count_wet, count_clear, dtype="float32")

    count_wet.attrs["nodata"] = nodata
    count_clear.attrs["nodata"] = nodata

    is_ok = count_some > 0
    count_wet = keep_good_only(count_wet, is_ok)
    count_clear = keep_good_only(count_clear, is_ok)

    return xr.Dataset(
        dict(
            count_wet=count_wet,
            count_clear=count_clear,
            frequency=frequency,
        )
    )
In [18]:
summary = reduce(wofl)

Wet counts¶

A count of every time a pixel was observed and wet.

In [19]:
summary.count_wet.plot(size=10)
plt.title("Lake Tempe – Wet observations counts for 2020");

Clear observations counts¶

A count of every time a pixel was observed (not obscured by terrain or clouds).

In [20]:
summary.count_clear.plot(size=10)
plt.title("Lake Tempe – Clear observations counts for 2020");

Wet frequencies¶

What fraction of the time was the pixel wet.

In [21]:
summary.frequency.plot(size=10)
plt.title("Lake Tempe – Wet frequency for 2020");
In [22]:
# Save frequencies to high-res image file
summary.frequency.plot(size=20)
plt.title("Lake Tempe – Wet frequency for 2020")
plt.savefig('/home/jovyan/tempe_wofs_2020.png', dpi=600);